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We present the application of a novel method of time-series analysis, the Hilbert-Huang Transform, 
to the search for gravitational waves. This algorithm is adaptive and does not impose a basis set on 
the data, and thus the time-frequency decomposition it provides is not limited by time-frequency 
uncertainty spreading. Because of its high time-frequency resolution it has important applications 
to both signal detection and instrumental characterization. Applications to the data analysis of the 
ground and space based gravitational wave detectors, LIGO and LISA, are described. 
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Introduction — A vital area of gravitational wave 
(GW) research is the development of data analysis al- 
gorithms that are able to detect non-linear, transient 
signals. The data analysis effort is especially relevant 
as the Laser Interferometer Gravitational Wave Obser- 
vatory (LIGO) [1, 2] undertakes a one-year science run 
at full design sensitivity. Also, the Laser Interferometer 
Space Antenna (LISA) [3] is engaged in the development 
of analysis algorithms to search for its candidate sources. 

The Hilbert Huang Transform (HHT) [4] presents a 
fundamentally new approach to the analysis of time se- 
ries data. Its essential feature is the use of an adap- 
tive time-frequency decomposition that does not impose 
a fixed basis set on the data, and therefore it is not lim- 
ited by the time- frequency uncertainty relation charac- 
teristic of Fourier or Wavelet analysis. This leads to a 
highly efficient tool for the investigation of transient and 
non-linear features. Applications of the HHT include ma- 
terials damage detection [5] and biomedical monitoring 
[6, 7]. Because General Relativity is an inherently non- 
linear theory, and because LIGO, LISA, and other GW 
detectors produce a great variety of non-linear and tran- 
sient signals, the HHT has the promise of being a power- 
ful new tool in the search for gravitational waves. This 
article describes the application of the HHT to GW data 
analysis. 

The HHT proceeds in two steps. First, the process of 
Empirical Mode Decomposition (EMD) reduces the time- 
series under analysis into components, known as Intrinsic 
Mode Functions (IMFs), thereby ’’sifting” or separating 
out the different frequency scales of the data. The IMFs 
when summed reproduce the original time series. The 
sifting is done adaptively, with no a priori structure im- 
posed on the data. 

The IMFs have a vertically symmetric and narrowband 
form that allow the second step of the HHT to be ap- 
plied: the Hilbert transform of each IMF. As explained 
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below, the Hilbert transform obtains the best fit of a si- 
nusoid to each IMF at every point in time, identifying an 
instantaneous frequency (IF), along with its associated 
instantaneous amplitude (IA). The IF and I A provide a 
time-frequency decomposition of the data that is highly 
effective at resolving non-linear and transient features. 

Instantaneous Frequency — The IF is generally ob- 
tained from the phase of a complex signal z(t) which is 
constructed by analytical continuation of the real signal 
x(i) onto the complex plane [7]. By definition, the ana- 
lytic signal is 

z{t) = x{t) 4- iy{t) (1) 

where y(t) is given by the Hilbert Transform: 

(Here P denotes the principal Cauchy value.) The am- 
plitude and phase of the analytic signal are defined in the 
usual manner: a(t) = \z(t)\ and 0(t) = ar gz(t) .The ana- 
lytic signal represents the time-series as a slowly varying 
amplitude envelope modulating a faster varying phase 
function [8]. The IF is then given by c o{t) = dO(t)/dt , 
while the I A is a (t) . We emphasize that the IF, a function 
of time, has a very different meaning from the Fourier 
frequency, which is constant across the data record being 
transformed. Indeed, as the IF is a continuous function, 
it may express a modulation of a base frequency over a 
small fraction of the base wave-cycle. Conditions on the 
form of x(t) that allow the application of this procedure 
are: 1) vertical symmetry with respect to the local zero 
mean, and 2) the same number of zero crossings and ex- 
trema [4]. 

A time-frequency analysis of a frequency modulated 
signal serves to illustrate the concept of instantaneous 
frequency. An important GW source for LISA is the 
inspiral of supermassive black holes (SMBH), The GW 
waveform from this source has a base frequency which 
chirps as the masses approach, and is also phase and 
amplitude modulated due to spin-orbit coupling [9]. The 
time evolution of the base frequency of the GW from the 



inspiral is given to first order by [10]: 
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where p = mim 2 /(mi+m 2 ) is the system reduced mass, 
M t is the system total mass, G is, Newton’s constant, c 
is the speed of light, and to is the coalescence time. 

Fig.l-(a) shows the waveform from an SMBH binary 
inspiral with masses m\ and m 2 of 10 6 and 10 5 so- 
lar mass, and spins $1 and S 2 of 0.95mf and 0.95m2. 
This waveform is used as input to the HHT and FFT- 
based spectrogram to perform the time-frequency anal- 
ysis, shown in Figs.l-(b), -(c), -(d). It is clear that the 
extraction of the IF and IA provides high time-frequency 
resolution of the waveform, showing detail of the chirp as 
predicted by Eq.(3) (dotted line in Fig.l-(b)), and also 
amplitude modulation due to the spin-orbit coupling. In 
contrast, the decomposition of the chirp through the FFT 
is seen to be limited by the time- frequency uncertainty re- 
lation: At A/ ~ I, where At and Af are the time interval 
and (Fourier) frequency bandwidth. The decoupling of 
the modulations of the IF and IA provided by the HHT 
allows a detailed look at the amplitude and frequency 
evolution of the waveform, and will aid significantly in 
determining the SMBH binary parameters. 

Empirical Mode Decomposition — The waveform of the 
previous example has an approximate symmetry about 
zero that satisfies the conditions listed above and allows 
meaningful application of the Hilbert Transform for ex- 
traction of the IF. To express an arbitrary time-series in 
this form the process of EMD is applied [4] . 

The EMD process consists of forming an envelope 
about the extrema of the data by cubic spline fitting, 
and then subtracting the average of the envelope from 
the data. The extrema of remainder are then fitted, and 
the process is repeated as many times as necessary to 
obtain a waveform that is symmetric about zero mean 
(within a predetermined tolerance). Once this has oc- 
curred, the waveform is labeled IMF1 and is subtracted 
from the original time series, removing the highest fre- 
quency content, and allowing the remainder to be sifted 
again to obtain the next IMF. This procedure is illus- 
trated in Fig. 2. 

The sifting identifies and removes components of the 
data first at the highest frequencies, then down in fre- 
quency to the lowest trends. For white noise the EMD 
acts as a dyadic filter, so that the central frequency as- 
sociated with IMFn is of order f s amp/2 n where / samp is 
the data sampling frequency [11]. A consequence of the 
sifting is that the IMFs are vertically symmetric about 
zero mean (although both the frequency and amplitude 
of an IMF vary in time.) The narrowband and symmet- 
ric form of the IMFs are essential to allow the Hilbert 
Transform to be applied in a meaningful way [4]. 

EMD and the Hilbert Transform: The HHT— To illus- 
trate the application of the full HHT to GW analysis, we 
look at the identification of a signal in white noise (Fig.3). 
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FIG. 1: Comparison of IF and FFT-based spectrogram of 
GW signal from SMBH inspiral with spin-orbit coupling, (a) 
GW waveform. (b),(c) IF and IA show sharp detail of fre- 
quency and amplitude modulation. The dotted line on (b) 
plots Eq.(3). (d) FFT spectrogram is blurred because of time- 
frequency spreading and spurious harmonics. 


As an example relevant to LIGO analysis, where strong 
signals are not expected, we consider here a weak signal 
relative to the noise. An issue for the HHT is that low 
signal to noise tends to degrade the IF and IA identifica- 
tion, as the signal extrema becomes distorted. A rough 
measure of this tendency was determined from simula- 
tions by injecting sine waves of varied signal to noise, 
frequency an d duratio n into white noise. This showed 
A/// ~ 2 \JStj 5msec / SN, where Af is the extracted 
frequency width of a periodic wave of frequency / and 
duration St in white noise, with signal to noise SN [12]. 
Thus in searching for signals with SN < 20 the HHT 
works best for St < 20msec. 

As an example we inject a 20 solar mass black hole 





FIG. 2: Illustration of Empirical Mode Decomposition (from 
Ref. [4]). (a) time series, (b) average of envelope formed by 
fitting extrema; (c) subtraction of average from time-series. 
Because the result is not vertically symmetric, the waveform 
will be re-sifted until an IMF is identified. 




FIG. 3: HHT of 20 solar mass black hole binary merger sig- 
nal in white noise, (a) merger signal, (b) signal in noise 
with SN=10. (c),(d),(e) IMFs 2,3,4. (f),(g) instantaneous 
frequency and power derived from the Hilbert Transform of 
IMF3 show the time-frequency-power structure of the merger. 


binary merger and ringdown signal [13] with SN=10 and 
a time duration of 5 msec into white noise at 16 kHz 
sampling rate. The merger signal is shown in Fig.3-(a) 
and the time-series of signal in noise is shown in Fig.3- 
(b). Figures 3-(c),-(d),~(e) show the sifted IMFs; in the 
3rd IMF the signal can be seen, largely separate from the 
noise. 

Figures 3-(f) and 3~(g) show the instantaneous fre- 
quency and power derived from the Hilbert Transform, 
of IMF3. The signal frequency and power are clearly 
visible in the Hilbert Transform of IMF3: the frequency 
shows a ramp during the merger and levels off to a con- 
stant during the ringdown, while the power shows the ex- 
pected rise and decay. This level of detail will aid signal 
identification, allowing comparison of signals from mul- 
tiple detectors. This may be contrasted with the Fourier 
Transform of the merger, which would give the power of 
only 1 point with a 200 Hz frequency spread for the en- 
tire 5 msec. Finally, we note that further noise reduction 
through averaging of the IF and IA over time is possible, 
as they are oversampled: the plots of Fig. 3 are sampled 
at the LIGO data rate of 16 kHz while the frequencies of 
interest for LIGO analysis are typically below 1 kHz. 

In analogy with an FFT-based search algorithm that 
looks for excess power in the Fourier power spectrum 
[14], the HHT can be used to search for excess power in 
the' time domain. For a given time record, the summed 
power is computed, and compared to the background 
level; thus clear evidence of a signal is seen in Fig. 3- 
(g), the power associated with IMF3. In comparison to 
an FFT-based search, the HHT is most effective for short- 
signals (< 20msec), where FFT analysis tends to lose sen- 
sitivity from time-frequency spreading. Thus the HHT is 
a useful tool, for example, in searching for black hole 
binary merger signals up to 100 solar masses. 

Analysis of Spectral Shapes — To facilitate the analysis 
of data spectral shapes, we use the HHT analogy of the 
power spectrum, called the ” marginal spectrum.” The 
marginal spectrum is defined as: 

M(cj) = ^ [ \ai(uj,t)\ 2 dt (4) 

i J 0 

where a*(a;,£) is the amplitude as a function of the in- 
stantaneous frequency at a given point in time and i is 
an IMF summation index. Thus the marginal spectrum 
is a measure of the total power present at a frequency to 
over the time interval T. 

Figure 4 shows the marginal spectrum of a time-series 
of length 1/16 second with the following Fourier power 
spectral shapes: fiat, proportional to frequency, and in- 
versely proportional to frequency. The marginal spectra 
are seen to have the expected shapes. Since the shape of 
the marginal spectrum does not depend on the integra- 
tion time, it may be used to examine the detector noise 
stationarity in fine detail. 

Application of HHT to GW detector characterization 
— The LIGO detector includes a large number of me- 
chanical resonators, servos, cavity optical resonances, 











FIG. 4: HHT marginal spectrum of l/16sec time-series with 
Fourier spectral shapes that are flat, proportional to fre- 
quency, and proportional to 1 /frequency. 
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FIG. 5: Time-frequency analysis of simulated instrumental 
artifact (Eq.(o)). a) Time-series of frequency modulated, 
40 msec transient in white noise, b) IF and IA from HHT: 
the modulation of the 60 Hz carrier from 20 Hz to 100 Hz 
is clearly visible, c) FFT-based spectrogram: time-frequency 
spreading blurs the modulation so that detail of its time struc- 
ture is lost. 


and light and sideband frequencies with differing Q’s and 
transient excitation levels [15], They can produce signals 
which can interact with each other in very transient and 
non-linear ways, for example upconverting noise at low 
frequencies into the signal band [16]. The HHT can be 
used to examine the detailed dynamics of their interac- 
tions in a way that is unavailable to Fourier decomposi- 
tion. 

We show the capability of the HHT to resolve the fol- 
lowing non-linear, transient waveform: 


y = cos [27r60t 4- sin (27r40t)] (5) 


which contains a 60 Hz carrier frequency, phase mod- 
ulated at 40 Hz. This could be produced by a laser 
frequency variation at 60 Hz coupled to seismic isola- 
tion stack motion at 40 Hz, for example. This waveform, 
which lasts only 40 msec, is shown adcled to white noise 
in Fig. 5- (a). 

After application of the HHT, the modulation of the 
carrier frequency from 20 Hz to 100 Hz is clearly visi- 
ble in Fig.5-(b), where the IA is shown as a function of 
color. In Fig.5-(b) the process of the upconversion of 
noise components at 40 and 60 Hz (outside the LIGO 
signal band), into components at 100 Hz (within the sig- 
nal band), is clearly revealed. Resolution at this level 
will aid in understanding transient experimental artifacts 
and distinguishing them from real signals. In contrast, 
Fig.5-(c) shows the FFT-based spectrogram, where time- 
frequency spreading blurs the modulation. In this figure 
the detail of the time structure of the modulation is lost. 

Summary — We have presented the application of the 
Hilbert-Huang Transform to the analysis of GW data. 
The high time-frequency resolution of the HHT, which 
is not limited by time-frequency spreading, will have im- 
portant application to both signal detection and detector 
characterization. 
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